# Replication File for "Combining List Experiment and Direct Question Estimates of Sensitive Behavior Prevalence"
# Peter M. Aronow, Alexander Coppock, Forrest W. Crawford, Donald P. Green
# Placebo Test II Power Calculations

rm(list=ls())

# Set your working directory to ACCG_list_replication_archive
setwd()

library(lattice)

#parameters for all figures
col.l <- colorRampPalette(c('white', 'black'))(100)
col.l[c(75:85)] <- c("red4","red3","red3","red2", "red2", "red","red2","red2","red3", "red3","red4")
depth.breaks <- do.breaks(c(0,1), 100)
pr.yes.treats <-  seq(0,1,by=.01)
pr.yes.controls <-  seq(0,1,by=.01)


# Figure A2: Power of Placebo Test II to Detect Violations of Treatment Independence

# Panel A2-1: N=100
powerz <- NULL
pr.yes.treatz <- NULL
pr.yes.controlz <- NULL

for(i in 1:length(pr.yes.treats)){
  for(j in 1:length(pr.yes.controls)){
    pr.yes.treat <- pr.yes.treats[i]
    pr.yes.control <- pr.yes.controls[j]
    power <- power.prop.test(n=50, p1=pr.yes.treat, p2=pr.yes.control)$power
    pr.yes.treatz <- c(pr.yes.treatz, pr.yes.treat)
    pr.yes.controlz <- c(pr.yes.controlz, pr.yes.control)
    powerz <- c(powerz, power)
  }
}

graph.frame.1 <- data.frame(x=pr.yes.treatz, y=pr.yes.controlz, z = powerz)

powerfig1 <- levelplot(z~x*y, graph.frame.1, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Pr(Y=1|Z=0)",
                       xlab = "Pr(Y=1|Z=1)",
                       main = "Power of Placebo Test II \n N= 100",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

# Panel A2-2: N=200
powerz <- NULL
pr.yes.treatz <- NULL
pr.yes.controlz <- NULL

for(i in 1:length(pr.yes.treats)){
  for(j in 1:length(pr.yes.controls)){
    pr.yes.treat <- pr.yes.treats[i]
    pr.yes.control <- pr.yes.controls[j]
    power <- power.prop.test(n=100, p1=pr.yes.treat, p2=pr.yes.control, strict=TRUE)$power
    pr.yes.treatz <- c(pr.yes.treatz, pr.yes.treat)
    pr.yes.controlz <- c(pr.yes.controlz, pr.yes.control)
    powerz <- c(powerz, power)
  }
}

graph.frame.2 <- data.frame(x=pr.yes.treatz, y=pr.yes.controlz, z = powerz)

powerfig2 <- levelplot(z~x*y, graph.frame.2, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Pr(Y=1|Z=0)",
                       xlab = "Pr(Y=1|Z=1)",
                       main = "Power of Placebo Test II \n N= 200",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

# Panel A2-3: N=500
powerz <- NULL
pr.yes.treatz <- NULL
pr.yes.controlz <- NULL

for(i in 1:length(pr.yes.treats)){
  for(j in 1:length(pr.yes.controls)){
    pr.yes.treat <- pr.yes.treats[i]
    pr.yes.control <- pr.yes.controls[j]
    power <- power.prop.test(n=250, p1=pr.yes.treat, p2=pr.yes.control, strict=TRUE)$power
    pr.yes.treatz <- c(pr.yes.treatz, pr.yes.treat)
    pr.yes.controlz <- c(pr.yes.controlz, pr.yes.control)
    powerz <- c(powerz, power)
  }
}

graph.frame.3 <- data.frame(x=pr.yes.treatz, y=pr.yes.controlz, z = powerz)

powerfig3 <- levelplot(z~x*y, graph.frame.3, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Pr(Y=1|Z=0)",
                       xlab = "Pr(Y=1|Z=1)",
                       main = "Power of Placebo Test II \n N= 500",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

# Panel A2-4: N=1000
powerz <- NULL
pr.yes.treatz <- NULL
pr.yes.controlz <- NULL

for(i in 1:length(pr.yes.treats)){
  for(j in 1:length(pr.yes.controls)){
    pr.yes.treat <- pr.yes.treats[i]
    pr.yes.control <- pr.yes.controls[j]
    power <- power.prop.test(n=500, p1=pr.yes.treat, p2=pr.yes.control, strict=TRUE)$power
    pr.yes.treatz <- c(pr.yes.treatz, pr.yes.treat)
    pr.yes.controlz <- c(pr.yes.controlz, pr.yes.control)
    powerz <- c(powerz, power)
  }
}

graph.frame.4 <- data.frame(x=pr.yes.treatz, y=pr.yes.controlz, z = powerz)

powerfig4 <- levelplot(z~x*y, graph.frame.4, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Pr(Y=1|Z=0)",
                       xlab = "Pr(Y=1|Z=1)",
                       main = "Power of Placebo Test II \n N= 1000",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

pdf("placebopowerIIfig1.pdf")
powerfig1
dev.off()

pdf("placebopowerIIfig2.pdf")
powerfig2
dev.off()

pdf("placebopowerIIfig3.pdf")
powerfig3
dev.off()

pdf("placebopowerIIfig4.pdf")
powerfig4
dev.off()

# Figure A2
pdf("powerIItiles.pdf", height=12, width=12)
print(powerfig1, split=c(1,1,2,2) , more=TRUE ) 
print(powerfig2, split=c(2,1,2,2) , more=TRUE ) 
print(powerfig3, split=c(1,2,2,2) , more=TRUE) 
print(powerfig4, split=c(2,2,2,2) ) 
dev.off()



